A mean field description of jamming in non— cohesive frictionless particulate systems 
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A theory for kinetic arrest in isotropic systems of repulsive, radially-interacting particles is pre- 
sented that predicts exponents for the scaling of various macroscopic quantities near the rigidity 
transition that are in agreement with simulations, including the non-trivial shear exponent. Both 
statics and dynamics are treated in a simplified, one-particle level description, and coupled via the 
assumption that kinetic arrest occurs on the boundary between mechanically stable and unstable 
regions of the static parameter diagram. This suggests the arrested states observed in simulations 
are at (or near) an elastic buckling transition. Some additional numerical evidence to confirm the 
scaling of microscopic quantities is also provided. 



I. INTRODUCTION 

Determining when and how a self-assembled partic- 
ulate system is capable of supporting macroscopic loads 
appears to be becoming established as a core paradigm in 
non-equilibrium condensed matter. For manydissipative 
athermal systems, such as granular media 0, y, C3- C3- 13- H3- 
Q , foams [a El E3 or emulsions [ill [l2| , the removal of a 
heat bath or other energy source results in a monotonic 
decrease of particle motion until a state of kinetic ar- 
rest is reached. The morphology, and hence mechanical 
response, of this static phase is determined by the dy- 
namics of its preparation Q , demonstrating the central 
importance of the initial dynamic phase and highlight- 
ing the non-equilibrium nature of the problem. It is this 
aspect that delineates self-assembled systems to those 
with a predefined geometry, such as the tensegrity struc- 
tures of great importance to engineering 114L Tl5| . or 
the disordered discrete and continuous models of rigidity 
percolation [H [H H H El El US- 

A generic feature of these distinct but related problems 
is the existence of a rigidity transition between states 
that have non-zero elastic moduli, and those that do 

not m mm m m m m hi. F 0r bond diluted 

lattices, the rigidity transition corresponds to a critical 
bond dilution, or equivalently a critical mean coordina- 
tion number z c , that is well approximated by the Maxwell 
constraint counting method [13j . This transition appears 
to bear some of the hallmarks of a continuous phase 
transition in equilibrium systems, including a diverging 
length scale, although other properties such as universal- 
ity remain unproven [2(j. For particulate systems, the 
picture is less clear. Central force models, such as fric- 
tionless elastic spheres in the limit of small deformation, 
appear to reach a similar transition 0, |29| , also with a 
diverging length scale |3(j , in the limit of infinite interac- 
tion stiffness under controlled pressure. If the volume is 
controlled, a critical volume fraction must be fine tuned. 
However, the introduction of friction, at least, appears to 
overshoot the critical z c by at least a few percent when 
in gravity 4, 5] , the exact amount apparently depending 
on the damping. 

Even when the transition can be approached arbitrar- 




FIG. 1: Combined dynamic flow and static stability diagrams 
in {z,u)) space, with z the mean coordination number and 
uj > a quantity proportional to the mean particle overlap. 
Solid (dashed) lines correspond to lines of constant volume 
V (pressure P), with the direction of minimising internal en- 
ergy U (enthalpy H) given by the arrows. Lines to the lower 
left corner correspond to smaller pressure or larger volume. 
The system becomes kinetically arrested at the boundary of 
the (shaded) stable region. The transition point at (z c ,cj) is 
shown as a solid disc (see text for discussion of z c ). 



ily closely, there remain crucial differences between dis- 
ordered lattices and particulate systems. In particulate 
systems, the morphology of the solid state depends on the 
dynamics of the immediate precursor to arrest, and may 
therefore include non-trivial, correlated structures. Most 
studies of lattice models assume morhpologies with no 
structure beyond the one-bond level. Furthermore parti- 
cles may become arrested in a state with a finite pressure 
P > 0, quite unlike the stressless configurations typically 
adopted in lattice models. Attempts have been made 
to circumvent some of the deficiencies of lattice models 
by postulating initial configuration generating algorithms 
that are hoped to mimic the dynamic self-assembly of 
particulate or atomic systems. Indeed one such scheme, 
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known as bootstrap percolation, has already been ap- 
plied to this problem |35j, and by infinite dimensional 
calculations claim to determine the exponents observed 
in simulations of central force systems. 

The purpose of this paper is to describe an analyt- 
ical scheme that, when applied to isotropic systems of 
particles interacting via central forces, derives the ex- 
ponents relating pressure, shear modulus and volume 
fraction to z — z c observed in simulations. There is 
no mapping to a known model, or assumption of any 
statistical mechanical analogy, as recently suggested by 
the elegant theory of Henkes et al. 36]. Instead, the 
properties of the static system are approximated accord- 
ing to a minimum assumption philosophy that proposes 
a series of intuitive approximations to locally close the 
equations. For the dynamics, a one-contact level clo- 
sure is employed, whereas for the statics the approxi- 
mation takes place somewhere between the one-contact 
and one-coordination shell level, in that it considers all 
contacts acting on a particle but ignores correlations be- 
tween them. The static aspect of this scheme has already 
been applied to mixed tensile-compressive systems |37^ . 

The central result of this paper is schematically repre- 
sented in Fig. ^ m which the statics and dynamics of the 
theory have been overlayed into a single plot. The axes 
are z and lo, where lu is proportional to the mean parti- 
cle overlap, so higher lu corresponds to increasingly com- 
pressive contacts. The rigidity transition lies at z = z c , 
lo = 0, which corresponds to the same unstressed transi- 
tion observed in e.g. lattice models. The anomalous low 
z c — 3 (which should be closer to 2c? = 6 in this d = 3 
dimension example) is explained later, but is essentially 
due to the reduction of degrees of freedom inherent in the 
approximation scheme. The shaded region corresponds 
to systems that are mechanically stable in that they gen- 
erate a positive restoring force in response to a localised, 
linear perturbation. The upper boundary of this stable 
region, which represents some form of buckling, is locally 
quadratic, and it is this shape that ultimately determines 
the exponents. 

This static picture is augmented by a global descrip- 
tion of the system dynamics from excited to arrested 
states. Again approximations are required, in this in- 
stance the kinetic energy is ignored and the system is 
assumed to evolve towards state with a lower internal 
energy (in the case of fixed volume) or enthalpy (in the 
case of fixed pressure), not unlike thermal systems |3S[ 
but without temperature/entropy terms. Some choice 
is required in defining volume in terms of z and lo, but 
for rapid quenches from highly excited states, a single- 
particle volume function should be valid, and will ro- 
bustly give the direction given in the figure. This drives 
the system towards the stable region. Kinetic arrest is 
assumed to take place when the boundary of the stable 
region has been reached; that is, the configurations ob- 
served in simulations lie on a line of buckling transitions. 
By suitably controlling pressure or volume, it is possible 
to bring the system to the same rigidity transition as in 



stressless systems, but along a different line to the lo = 
systems considered in disordered lattices. 



II. STATICS: THE MEAN MODE 
APPROXIMATION 

A variety of approximate analytical treatments for pre- 
dicting the mechanical response of athermal materials 
have been devised; just two will be mentioned here by 
way of comparison. Perhaps the most straightforward is 
to assume the induced deformation field is affine, so the 
interparticle displacements are just scaled-down versions 
of the macroscopic field. This has been applied to fric- 
tional granular media [3fll | , but is incapable of predicting 
a rigidity transition at finite density or volume fraction. 
Even side-stepping this deficiency by directly comparing 
elastic moduli to pressure does not resolve its failings [s3 j 
indicating that modes near the transition are inherently 
non-affine. Enhanced versions have been proposed in 
which additional degrees of freedom are introduced and 
determined by variational principles |4fl Ell ] , but only at 
the cost of significant additional complexity and still not 
suitable for studying the transition. 

Another approach, extensively applied to unstressed 
systems of predefined morphology, is known as EMA 
or the effective medium approximation (although this 
term is sometimes also applied to the affine approxima- 
tion Here the disordered system is replaced by a ho- 
mogenous analogue with effective interaction parameters. 
The procedure can be crudely summarised as follows: a 
single disordered element is inserted into the homogenous 
bulk and the response to this isolated defect determined 
using the Green's function. The effective parameters are 
determined by demanding that the mean response aver- 
aged over all possible states of the inserted element is 
zero. Although this has been successf ully ap plie d to a 
range of continuous and lattice systems [iH HH l42l Eflj . it 
has a fatal shortcoming when applied to particulate sys- 
tems: the calculation cannot proceed without first identi- 
fying an analogous homogeneous material with a known 
Green's function. It is difficult to postulate such an ana- 
logue for finite-size particles, except in certain special 
cases such as quasi-ordered systems Q. Furthermore 
we have already argued that prestress is likely to be im- 
portant, and even for lattice spring networks the Green's 
function with prestresses is not known exactly |24j . For 
granular media the problem is particularly acute, as the 
Green's response has been the subject of substantial de- 
bate in recent years (see e.g. 0,00 and refs therein). 
The approximation scheme detailed below circumvents 
these issues, as we describe below. 

For this first exposition of the theory, the non-damping 
part of the interparticle interaction is assumed to be a 
purely radial pair potential, generating central forces act- 
ing along lines connecting particle centres. This choice 
is made to reduce the parameter space to a manageable 
size and to allow a systematic and lucid unfolding of the 
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resulting phase structure of the system. Nonetheless this 
simplification should hold approximately true for emul- 
sions or wet foams near the transition, when all particles 
are only slightly deformed, although it will clearly not ap- 
ply to strongly deformed particles (or dry foams 8] ) . For 
real granular media there will also arise transverse surface 
forces mediated by non-zero friction at the bead-bead 
interface; purely central forces correspond to frictionless 
spherical beads, which only truly exist in numerical sim- 
ulations. 



A. Determining the mechanical stability 

Given central force interactions, there is no need to 
track particle orientations and the static system is fully 
specified by the (i-dimensional position vectors x' 3 of all 
particles (3. The force on (3 due to 7 is denoted by f 7 ^, 

fr/s _ /( r 7 ' 3 )n 7/3 (1) 

where r 7/3 =| x! 3 — x 7 | is the distance between particle 
centres and n 7 ^ = (x* 3 — x 7 )/r 7/3 is the unit vector from 
7 to (3. In this context the scalar central force f(r) is 
usually taken to be of the form 

/M^IK 1 ^)" 1 r<r ° (2) 
[ : otherwise 

with a = 1 (truncated Hookcan) or a = 3/2 (Hertzian), 
and ro is the sum of the two particle radii. The prefac- 
tor /i > is typically treated as a particle-independent 
parameter, although strictly speaking it is a function of 
the radii of curvature at the point of contact for Hertzian 
interactions [33l l45| . Note that with this sign convention, 
positive / corresponds to compressive forces and negative 
/ to tensile ones. 

Suppose we are given an initial configuration {x^} that 
is static, i.e. the vector sum of all contact forces on 
each particle vanishes. To determine its linear stability 
apply an arbitrarily small external force (5f cxt onto the 
particle lying nearest some arbitrary point in space - 
call this particle a (not to be confused with the force 
law exponent). If the system is stable, it will move to a 
nearby static configuration in which all particles (3 have 
been displaced to x /3 + 5x /3 with |<5x| <C r . Force balance 
must again be obeyed, i.e. the changes in contact forces 
on (3 sum to zero for (3 ^ a, and to — <5f cxt for particle a. 

The response <5x a for a particular configuration {x* 3 }, 
even if tractable, would be of no practical interest and 
we must instead ensemble average, keeping fixed a set 
of parameters that dominate the mechanical response of 
the system. A wealth of data has shown the mean coor- 
dination number z to be a crucial factor in determining 
stability, and Alexander nas highlighted the impor- 
tance of prestresses, so we also assume that both z and 
some measure of the initial contact force distribution is 



kept fixed. Once the ensemble (denoted by the angled 
brackets (. . .) below) has been suitably defined, the re- 
quirement of force balance on the perturbed bead a can 
be written 




where the sum is over all (3 interacting with a. The 
change in contact force Sf a @ can be related to the particle 
displacements Sx a and Sx 13 by 

5ff = Al?(6 x e-5x«) (4) 

with summation over Roman indices only. The d x d 
matrix A a @ is defined by ^(| 



Aij = —^—(Sij - hihj) + f'{r)n. l h j (5) 

assuming /(r) is continuous with a finite first derivative 
f'(r) 7^ over all r of interest. Here and below the 
suffices a, (3 are dropped whenever the meaning is clear. 



B. Derivation of the MM A 

So far this is exact but intractable. Now we approx- 
imate. Consider that, for an isotropic system, the per- 
turbed bead must move parallel to the external force af- 
ter averaging, (5x a ) = A<5f cxt with an unknown compli- 
ance A. The philosophy of the MMA is to impose this 
form before averaging, i.e. inside the brackets in 
In this way the dependency of Sx a on the entire initial 
configuration {x 7 } is subsumed into the single scalar pa- 
rameter A. This is clearly a significant saving in terms of 
complexity, although it disallows the transverse motion 
of the particle and hence reduces the degrees of freedom; 
the consequence of this on the location of the rigidity 
transition will be discussed later. 

The logical continuation of this approach is to similarly 
replace the Sx 13 by (Sx 13 ); however, this averaged form 
cannot be determined by symmetry considerations alone. 
Instead we assume here that the change in the contact 
force with a can be treated as an external force on (3, so 
that 6x? = XSi af3 with the same A as before. Intuitively, 
this corresponds to the statement that the displacement 
of (3 is dominated by the change in contact force with a, 
which, for a monotonically decaying force field extending 
outwards from a, should at least not be embarrassingly 
wrong. These two approximations taken together allows 
each contact force 5f a @ to be uniquely determined from 
(5f cxt , as found by inserting 8x a — A<5f cxt and Sx 13 = 
\5f a/3 into (0J and (JSJ and inverting, 
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Sff = S?P5f, 



s: 



a/3 



i] "J] ' 

l + (X\f'(r^)\Y 



-1 . a/3 ~a/3 



1 - 



a/3N 



r a/3 



fa 



~a(3 ~a(3 



(6) 



Thus each (5f Q/3 is independent of the others. Note that 
the unphysical singularity at A//r = 1 is avoided by the 
stability equation below. 

It is apparent from © that the MMA has reduced the 
global problem to a local one, in which the response of 
each contact Si a/3 depends only on the interparticle sepa- 
ration r a/3 (through which /(r) and f'(r) are found), the 
unit vector h a P , and the coordination number for this 
bead z Q , implicit in the summation 10. Before the av- 
eraging can be completed, it is necessary to specify how 
these quantities vary. Here we will deliberately take sim- 
ple forms to facilitate transparent interpretation of the 
results. Firstly, z a is taken to be independent of the con- 
tact forces and orientations, allowing the ^"-averaging 
to be performed and the the force balance equation © 
rewritten as 



s 



a/3 



a/3 



= 



(7) 



where z is the mean coordination number, and the aver- 
aging is now over n a P and r a/3 . Since <5f cxt is arbitrary, 
the quantity inside the brackets in Q must vanish. 

We further assume that the bond orientations are in- 
dependent of the contact forces. The n"' 3 are taken to 
be independent, identically distributed random variables, 
uniformly distributed over the d — 1 dimensional unit hy- 
persphere. This neglects any correlations in the topology 
of the network. It also ignores excluded volume effects, as 
it allows the same particle to have bonds arbitrarily close 
together (so the contacting particles would significantly 
overlap); this should not be crucial near the transition, 
which is determined by stability requirements rather than 
excluded volume, but will be relevant at higher z. Per- 
forming the average gives 



i[i -(A/TT 1 



= 1 



(8) 



where the identity {fiifij) = 4<Sy has been used. 

It remains to specify how the contact forces are dis- 
tributed. In principle only configurations consistent with 
force balance should be allowed, but this complication 
becomes redundant given the approximations leading 
to JJjJ, which allows the response from each contact to 
be calculated independently. A natural choice is then 



to assume that each contact force f(r), or equivalently 
each interparticle separation r, are identically and inde- 
pendently distributed according to some given distribu- 
tion. Clearly this neglects any correlations in the initial 
force network, but force balance is ensured on average 
by virtue of the uniform distribution of contact angles 
already employed. Some calculations for general force 
distributions will be described later. For now, the sim- 
plest choice possible is made, namely a delta-function 
distribution corresponding to a monodisperse separation 
r, force f(r) and gradient f'(r) ^ for each contact. It 
is then possible to insert @J into and integrate; the 
result is finally 



d[--l)=(d- 



1) 



1 



A/(r 



- 1 



1 + A | /' | 



(9) 



This is a quadratic equation in A, which is easily solved. 

Although A is in principle a measurable quantity, it is 
more useful to specify results in terms of the shear mod- 
ulus G or bulk modulus K. These can be related to the 
compliance A by specifying some suitable closed surface 
S and applying an external force to each particle it cuts, 
i.e. forces parallel to S for G and normal for K. Given 
the displacement of each particle is oc A -1 according to 
the MMA, it is straightforward to see that 



G ~ B G \- 



K ~ B K \~ 



(10) 



where the prefactors B have dimension (length) 2_<i . The 
relevant length scales are the characteristic length of the 
enclosing surface L ~ d ~^/~S and the particle radius r /2, 
but the local closure of the MMA equations means we 
are unable to determine their weighting in Bq and Bk- 
However, if A diverges with L, as is trivially true for a 
d = 1 system with fixed boundaries, then B ~ LrQ~ d to 
ensure finite moduli for arbitrarily large systems. Some 
form of divergence of A with L is also expected from linear 
continuum elasticity ^]|- This subtlety is sidestepped 
below, where L is assumed to be fixed and finite. 



C. Stability regimes 

Ignoring the trivial d — 1 case, dimensionality only en- 
ters via the prefactors and so all d > 2 will be discussed 
together. The solutions to © can be conveniently ex- 
pressed in terms of the two scalars z and u>, where 



f(r)/r 
I f'(r) | 
1/ _r 
a \ r 



(11) 
(12) 



is a dimensionless measure of the prestress in the system. 
The second form <|12[) holds for the particulate potentials 
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FIG. 2: Stable regimes of (z, u) space in the MMA model, 
for d = 3. The black disc at (d, 0) is the (unstressed) rigidity 
percolation transition. All points z < d, lo — are unstable. 
For lo < 0, corresponding to tensile bonds, the system is al- 
ways stable (light grey region). For lo > 0, corresponding to 
compression, only the dark grey shaded region is stable. 

PJl in the limit r — ► Tq, i.e. close to the rigidity tran- 
sition, which is the regime of interest here. A schematic 
description of the predictions of the MMA are given in 
Fig. [21 Some important features are now discussed. 

lo = 0: This is the unstressed case in which all con- 
tact forces are initially zero, nullifying use of the par- 
ticulate potentials J5J, which have no linear response at 
r = ro, but still attainable for non-truncated Hookean 
springs. The equation for A gives the single solution 
A|/'| = (z/d-l)-\ oyG,K ~ (z/z c -l)f with the tran- 
sition point z c — d and an exponent / = 1. The effective 
medium theory for diluted spring lattices also predicts 
/ = 1 , but at the higher transition point Zr = 2d in ac- 
cord with the Maxwell counting estimate |43 ■ The tran- 
sition value found here, z c = d, seems anomalous until 
one recalls that the basic assumptions of the MMA re- 
strict the motion of the particles to mean forms, thereby 
reducing their degrees of freedom and hence lowering z c . 
Despite this, the MMA still predicts a finite transition, 
and can therefore can be used qualitatively. Any unease 
over the actual value could be lessened by referring to it 
as an effective coordination number z if desired. 

lo < 0: When all of the bonds are tensile, there is always 
one real, positive solution of A extending from arbitrar- 
ily large z down to a lower value z m j n = 1. Again this 
value is too small; z m i n = 2 is more probable, i.e. infinite 
chains of particles spanning the system. As lo — > with 
z > z c fixed, the single root of A continuously approaches 
the unstressed solution given above. Repeating this pro- 
cedure for z < z c , however, reveals that A diverges as 
\lo\~ 1 and hence G,K ~ \lo\ vanishes continuously as the 



unstressed axis is reached. Thus just below the lo = 
line, the elastic moduli are very small and the system is 
inherently weak, becoming weaker as z decreases. This 
may explain why the few attempts to survey this region 
in disordered lattices [24l |25| have observed a rapid but 
gradual crossover of the transition from z c to z m - ln : nu- 
merical noise and/or arithmetic precision may incorrectly 
attribute zero values to small but finite moduli. 

lo > 0: For compressed bonds, the {z,lo) plane is parti- 
tioned into a stable region with two distinct real, positive 
roots, and an unstable region for which both roots are 
either complex or negative. The boundary between the 
stable and unstable regions is quadratic near z c , 



" bdy " Ad^d-iy z > Zc - (13) 

Both roots of A coincide on the boundary, 

^ z-z c r { 4dHd-l) \ a 
Abdy ~ 2d(z-l)j[a(z-z c ) 2 ) (M) 

and hence Gbdy,-Kbdy ~ Aj^ d ~ [z — z c ) 2 " -1 - Starting 
from the stable regime and decreasing lo to zero, one of 
the roots diverges as lo~ x while the other continuously ap- 
proaches the unstressed solution, crossing over to become 
the single root in the tensile regime (where the other root 
becomes negative). 

The manner in which the compressive system becomes 
unstable is noteworthy. On the boundary, ro — r ~ 
(z - z c ) 2 , / ~ (z- z c ) 2a and /' ~ (z - Zc) 2 ^" 1 ), which 
according the means that the force transfer is predom- 
inantly longitudinal. As already noted by Alexander |47| . 
in such cases the change in energy will be positive, from 
which he infers the system should be stable. However, 
there are other ways of buckling. An established alterna- 
tive is a bifurcation to a different class of solution [48| : 
we might also speculate that the energy landscape may 
exhibit discontinuities in the limit of infinite system size, 
allowing some form of catastrophic buckling. In fact, the 
buckling as envisaged by Alexander, which corresponds 
to A < here, does arise within the MMA, but only for 
z < z c and small lo > 0. The upper boundary in Fig. [21 
rather corresponds to when A becomes complex. 

III. DYNAMICS: ENERGY MINIMISATION 

A system not in a shaded region in Fig. [21 will desta- 
bilise under any non-zero noise, evolving its contact net- 
work according to the dynamical particle interactions and 
hence allowing z and lo to vary. It can only come to rest 
in a mechanically stable region. Indeed for sufficiently 
damped interactions, as assumed here, kinetic arrest can 
be identified with the point at which the system first 
touches a stable region. Strong damping also means that 
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the kinetic energy is always small, so that the system will 
evolve to minimise some suitably defined energy poten- 
tial. For constant volume V, this potential is the internal 
energy U , here just the total potential energy stored in 
the interparticle bonds. For controlled pressure, the cor- 
responding potential is the enthalpy H — U + PV 38] . 

A crucial problem in integrating the fixed P or V dy- 
namics and the (z,u>) stability diagram is writing down 
expressions relating V to z and lu. This is likely to be 
a subtle issue; under gentle shaking, the particles may 
form spatially extended structures that would necessitate 
many-particle variables to calculate V Q, IHJ] , which 
is clearly beyond the one-particle closure of the MMA 
equations. For now we ignore such potential pitfalls, and 
instead assume the following, one-particle description, in 
the expectation that it will hold in the initial dynamic 
phase from a highly excited initial state. Simply assume 
that V is a decreasing function of both z and u>, as might 
be expected for uniform, global changes of these vari- 
ables. The precise choice of V(z, w) should incorporate 
the large changes in z that are possible for small changes 
in r when the particles are barely touching, i.e. when 
lu <C 1. This can be written as 



Dlu" 



(15) 



where the unknown exponent b is assumed here to obey 
b > 1 (b < 1 alters the scaling behaviour of V with z 
described later, but not of P, G or K). The dimensionless 
constant D > is some material-dependent parameter. 

It is straightforward to derive the internal energy U by 
employing the same approximations as used to perform 
the integration of the MMA equations, namely a constant 
u — —(1 — r/ro) and isotropic bond orientations n, 



U 



Nz tq/x 



T (^)« +1 



(16) 



which is the total number of contacts Nz/2 multiplied 
by the bond potential. Similarly, the isotropic pressure 
PSij is the sum of the r^fj for each bond, divided by the 
volume V, or (for r w ro) 



(17) 



where the identity {fiifij) = has been used. 

Performing the minimisation for both cases reveals 
broadly the same behaviour; the system will evolve in 
the direction of increasing z and decreasing to, as schema- 
tised in Fig. ^ and already discussed in the introduction. 
For example, for fixed V, the extremum of U (which we 
assume is the minimum) is found by solving dU = si- 
multaneously with dV = 0, where the latter gives the 
constraint of constant volume. This can be rearranged 
to give U, w V, z — U, z V,u) and hence from ljTU|l . 



zto a V, z = 



,a+l 



a 



1 



(18) 



This admits the single solution lu — in the small-w 
regime of interest here, so U is minimised when all parti- 
cles are at the limits of their interaction potentials (or 
at their natural lengths for Hookean springs). Given 
that the system is constrained to move on lines of con- 
stant V (I15f) . the minimum corresponds to a divergent z, 
clearly unobtainable in a real system. Excluded volume 
and ordering effects must be incorporated into the theory 
before any large-z treatment can be attempted. Repeat- 
ing the enthalpy minimisation at fixed P gives essentially 
the same behaviour. 



A. Kinetic arrest and scaling behaviour 

Given minimisation drives the system in the direction 
of small lo and large z, they will enter the mechanically 
stable region of the (z,lu) diagram Fig. ^ somewhere 
along the stability boundary. For overdamped dynam- 
ics, we assume it also stops there, allowing the scaling 
behaviour of various quantities with z — z c to be de- 
termined. According to (12} , the microscopic variable 
u> r~j (z — z c ) 2 , which is confirmed by d — 3 numerical 
simulation of Hertzian spheres shown in Fig. [3] The re- 
lation between the elastic modulus G and z — z c on the 
boundary has already been given (|14H . 



g ~ (z - Zc y 



(19) 



i + 



with c = 2a — 1. Note that this diverges as a , 
signifying the breakdown of linear response in this ad- 
mittedly atypical class of pair potentials. At the transi- 
tion point (z,lu) = (z c ,0), the pressure P is zero with a 
finite volume Vo , so PV « PVq to first order and scaling 
relations involving P can be found independently of the 
choice of V (fT5|l . From fT7|l . 



P~ {z-ZcY 



(20) 



with e — 2a. Finally, for relations involving V , and hence 
the volume fraction <fi ~ V^ 1 , we find 



v -v~(z- z c y 



(21) 



with 5 = 2 when b > 1 (g = 2b for b < 1). These results 
are summarised in Table [lj where they are compared to 
simulation results on various central force systems. 

It should be noted that the postulated forms for Vo — V 
and P already allow the bulk modulus K to be deter- 
mined without recourse to calculating the mechanical re- 
sponse, and as expected (given the one-bond expressions 
for P and V), the predicted exponent is that of an affine 
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FIG. 3: Confirmation of the predicted scaling (uj) ~ (z — z c ) 2 
for d = 3 Hertzian spheres relaxed at constant volume. The 
straight line is a fit to w a (z - z c ) h with h = 1.92(6) and 
z c = 5.86(4), where figures in brackets denote single stan- 
dard errors on the last digit and z c < 2d — 6 since (rigidly) 
disconnected rattler particles were included in the average. 
The simulations were performed on an Af = 1000 monodis- 
perse particle system using the conjugate gradient minimi- 
sation routine in the COGNAC module of the open source 
numerical suite OCTA ifl. 



deformation, K ~ P/(V - V) ~ (z - z c ) 2a ~ 2 , different 
to that of G given above. These different exponents have 
also been seen in simulations [2jJ, although it counters 
the consensus of lattice models, where all elastic mod- 
uli obey the same scaling behaviour is in mm. it 

also disagrees with the MMA exponents presented ear- 
lier l|10l) . It is not yet clear if this is a real physical 
phenomenon or a coincidence of anomalous simulation 
results and naive choices for P and V in this theory. In 
any case, we are forced to conclude that the MMA pre- 
dictions only apply to non-affme deformations. 

The similarity of the MMA exponent for G and that 
of simulations is noteworthy, as the theory relied on a 
local closure of equations and hence cannot incorporate 
the long wavelength modes observed in simulations [30| . 
Persisting with the continuous phase transition analogy, 
it is possible that the scaling regime for non-trivial expo- 
nents is incredibly narrow; this is not without precedent, 
as ionic fluids exhibit mean field exponents except very 
near the critical temperature [5 4. l55| . Alternatively the 
upper critical dimension may simply be 2. It should be 
noted that most (but not all |29f) simulations maintain 
a fixed system size, which will also give mean field ex- 
ponents when the correlation length exceeds the system 
size. 

A brief survey of other simulations not shown in the 
tabic indicate what new features will alter the scaling be- 
haviour and hence represent relevant perturbations. Fric- 
tion is an obvious candidate, and indeed recent simula- 
tions of Zhang et al. |12| demonstrate that infinite friction 
will alter the exponents. They also show that finite fric- 



tion introduces history dependency, suggesting a more 
advanced description of the dynamical phase will be re- 
quired for a full theory. Also, the molecular dynamics 
simulations of Kasahara and Nakanishi 0, Q appear to 
find exponents consistent with simple rationals that are 
nonetheless different to those predicted here. This may 
be because their system has gravity, and hence the con- 
tact network is anisotropic and the force balance equa- 
tions couple directly to an external field, either of which 
may be a relevant perturbation. 



B. Distributed contact forces 

Of all the enhancements to the MMA theory that could 
be incorporated, perhaps the most pressing is to relax the 
assumption that every contact force /, or equivalcntly 
every particle overlap 8 = r$ — r, is the same. Simula- 
tions have demonstrated that these distributions are in 
fact continuousl y di stributed right down to zero forces or 
overlaps [5(J, loTl l52l | . Below we present some calculations 
that probe the effect of polydisperse 8 within the MMA 
framework. Our conclusion is that it in fact makes very 
little difference to the overall behaviour of the system, 
and does not change the exponents already quoted. We 
are also able to confirm the scalin g of the overlap distri- 
bution as observed in simulations [29|. 

The simplest way to incorporate distributed overlaps 
is to retain the approximations leading to (JSJ, namely 
independent isotropic bond orientations n and a local 
coordination number independent of the contact forces, 
and perform the integral over a known, fixed distribu- 
tion P(5). This can be performed explicitly for a con- 
venient choice of parameters, for instance Hookean in- 
teractions a = 1 with a uniform overlap distribution 
P(S) — 1/S for < 6 < So- This then produces an 
equation for A that collapses into the form already stud- 
ied lO in the limit So — > 0, with the separation r replaced 
by a mean overlap ro — Sq/2. Given the volume function 
(|15|l still holds, we see that 5 being distributed uniformly 
down to 5 = does not significantly alter the statics or 
dynamics of the system in this case, but merely modifies 
the prefactors. 

More general distributions can be considered by as- 
suming that the behaviour observed for the monodisperse 
case still broadly applies. Specifically, we assume that a 
unique stressless rigidity transition exists at some point 
z = z c , and a boundary between stable and unstable 
compressive regimes extends continuously from this point 
into the region z > z c , similar to Fig. |2 Furthermore, 
we assume that the system becomes kinetically arrested 
on this boundary under energy minimisation. Therefore 
the averaged compliance A is expected to scale with the 
distance from the transition e = (z — z c )/z c > as 



A = A e~" (22) 
for small e, with v an unknown positive exponent. Then 
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we make the following scaling ansatz for the distribution 
of overlaps P(S), 



P(6) = E-^q(s-^S) 



(23) 



where q(x) is a fixed distribution. I)23|) states that the 
distribution of overlaps will uniformly contract as the 
transition is approached, with a width s that vanishes as 
s~£ 7 with 7 > 0. 

Inserting l)22|l and 1)23(1 into JHJ allows the integration 
to be performed. Different results occur for different com- 
binations of exponents, but only the form of solution for 
a"/ > v and 7(0; — 1) < v admits a stressless rigidity 
transition. The equation for Ao in this case is 



l- d - 

z 



(d- 



D^ 0fJ, c~ v+c 



ro 

_i/-7(q-1) / 1-a 



lq{x) 



\ /q(x) 



(24) 



where the angled brackets here denote averaging over the 
fixed distribution q(x). I|24|l admits a solution e = at 
z = z c = d given the exponent inequalities just quoted. 
It is similar to the monodisperse expression @ with / 
replaced by {x a ) q ^ x y and (x 1_a ) . , the quantity related 

to/'. 

Stable solutions of l|24|) correspond to Ao real and pos- 
itive, and so the boundary between the stable and un- 
stable regions can be found using the familiar quadratic 
equation formulae. For consistency with the earlier as- 
sumption that the compressed stability boundary for 
(x a ) q , x \ > is continuously connected to the point 
z — z c , £ — 0, we find it is necessary to impose 7 = 2. 
This means that the width of the overlap distribution 
scales as s ~ (z — z c ) 2 near the transition. The simula- 
tions of O'Hern et al. have shown that s ~ ((f)— C ) A with 
A close to 1 for Hookean interactions a = 1 in d = 3 [2^| . 
According to Tabled this corresponds to s ~ (z— z c ) 2A , 
in agreement with the prediction 7 = 2A = 2 found here. 
A second consistency check is that Ao = 0(1), which al- 
lows the second exponent to be fixed, v = 2a — 1. Note 
that both of these exponents are equal to their counter- 
parts in the monodisperse-i5 case, i.e. I|13l) for 7 and 
(|14|l for v. Although this analysis in non-rigourous, it 
strongly suggests that polydisperse contacts does not al- 
ter the scaling picture presented earlier. 



IV. DISCUSSION 

A central feature of this class of problem is the in- 
trinsic interweaving between the dissipative dynamics 
as the system cools, with the mechanical response of 
the arrested state. In this paper, a minimal coupling 
has been presented, namely that the dynamics proceeds 
in an independent -particle manner until a mechanically 



TABLE I: Table of the scaling relations between Az = 2 — z c , 
G, P and AV — Vo~ V ~ (4>—(f>o) (with <j>o the critical volume 
fraction) as predicted by the MMA theory, where a is the force 
law exponent @, d dimension. For comparison, results from 
simulations of central force systems and the trivial (affine) 
predictions are also shown. 



Model 


G ~ Az c P ~ Az e AV ~ Az 9 


MMA 
a > 0, d > 2 


2a - 1 2a 2 


Affine (see [29]) 


2a - 2 2a 


Wet foam [9] 
a= 1, d = 2 


«l a 2 ±0.4 2 ±0.4 


O'Hern et al. [291 

a = 1, d = 2,3 
a = 3/2, d = 2,3 


1.01 ±0.1° 2.1 ±0.2 2.04 ±0.1 
2.08 ±0.1 a 3.15 ± 0.3 2.08 ±0.1 


Zhang et al. [12] 
a = 1.28, d = 3 


w 2.45 « 1.96 


Makse et al. [3J 6 
a = 3/2, d = 3 


3.3 ±0.5 2.1 ±0.6 



"Result for shear modulus shown. 
''Only frictionless data shown. 



stable region is reached, Fig. ^ This is highly sim- 
plified, and a more elaborate theory is desirable, per- 
haps alon g th e lines of the bootstrap percolation model 
approach |35|. Intuitively, we expect that, during the 
dynamic phase, transient overconstrained, rigid clusters 
will become stressed by interparticle collisions, and re- 
lieve this stress by becoming non-rigid, i.e. expanding 
into neighbouring, underconstrained regions. Kinetic ar- 
rest occurs when a spanning rigid cluster forms. Note 
that this argument suggests a dynamic homogenisation 
process, perhaps explaining the apparent appearance of 
mean field exponents in d — 2 and 3 simulations in ta- 
bled 

This first application of the mean mode approxima- 
tion has been applied to arguably the simplest par- 
ticulate problem, namely repulsive central forces in an 
isotropic system. The simplicity of its results suggests 
that additional features could be included while remain- 
ing tractable. For instance, friction and gravity would be 
needed before any sensible comparison with real granu- 
lar media could be made. A problem that may emerge is 
closing the equations; the averaged force balance equa- 
tion J3J| only gives one scalar equation, which is why the 
proposed displacement modes were parameterised by sin- 
gle scalar A. If a future application had too few equations, 
one possible approach would be to assume the response 
will minimise the increase in elastic energy, converting it 
to a minimisation problem with any known equations as 
constraints. In principle, displacement modes with any 
number of unknown parameters could be introduced by 
this approach. 
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